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Abstract: Bacteria in desert soil have unique phylogeny and important ecological functions, and their 
responses to changes in precipitation need further attention. However, relevant studies have mainly 
focused on the surface soil, and studies on the responses of bacteria at different soil depths to variations 
in precipitation are rare. Thus, we used 16S rDNA high-throughput sequencing to investigate the changes 
in soil bacterial distribution along a mean annual precipitation gradient (50-150 mm) in the Alxa Desert, 
China, and compared the variation characteristics in the surface soil layer (0—10 cm) and subsurface soi 
layer (10-20 cm). Results showed that soil bacterial communities significantly changed along the 
precipitation gradient in both soil layers. However, the subsurface soil layer could support bacteria 
communities with higher diversity and closer internal relationships but more internal competition than the 
surface soil layer. Additionally, compared with the surface soil layer, variations in diversity and 
co-occurrence patterns in the subsurface soil layer were more in line with the changes in the mean annua 
precipitation, while bacterial community structure was less variable in the subsurface soil layer. Compared 
with the mean annual precipitation, soil moisture had little influence on the structure and diversity of soil 
bacterial community but had a high correlation with intercommunity connectivity. Therefore, soi 
moisture might play a complex role in mediating environmental conditions and soil bacterial community 
characteristics. Due to the different responses of surface and subsurface soil bacteria to the changes in 
precipitation, it is necessary to distinguish different soil layers when predicting the trends in desert soi 


bacterial conditions associated with precipitation, and prediction of subsurface soil bacteria may be more 
accurate. 
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1 Introduction 


Compared with humid climate zones, arid regions are more sensitive and vulnerable to climate 
change (Yi et al., 2014). Especially for middle and low latitude deserts, changes in precipitation 
may be the environmental factor with the greatest impact on ecosystem dynamics (Loarie et al., 
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2010; Yi et al., 2014). Furthermore, changes in water conditions induced by climate change will 
affect the structure and function of soil microbial community (Azua-bustos et al., 2012; Tomiolo 
et al., 2015), and such changes can affect the release or storage of soil carbon (Cregger et al., 
2014). Due to the severe environmental constraints on higher plant and animal life in the desert, 
soil microbial communities are considered to play a more critical ecological role (Pointing and 
Belnap, 2012) and have more potential to maintain the stability of soil ecosystem in the desert 
(Cruz-Martinez et al., 2009). Moreover, soil microbes are affected by both deterministic and 
nondeterministic factors (Zhou et al., 2004, Johnson et al., 2017), and the response degrees of soil 
microbe to changes in moisture differ (Manzoni et al., 2012). Therefore, before predicting how 
desert ecosystems will change in the context of future climate change, it is necessary to continue 
to deepen the understanding of the impact of precipitation changes on desert soil microorganisms. 

In view of the phylogenetic and functional peculiarities of desert soil microorganisms (Fierer et 
al., 2012), many studies have addressed the impacts of changes in precipitation or water 
conditions on these organisms (Clark et al., 2009; Angel et al., 2010; Bachar et al., 2010; Bachar 
et al., 2012; Barnard et al., 2013; Huang et al., 2015; Maestre et al., 2015; Tripathi et al., 2017; 
Zhang et al., 2018). Moreover, since soil microbes and organic carbon inputs are generally 
concentrated on the soil surface (Zvyagintsev et al., 1994), most of these studies concentrated on 
the soil surface or did not distinguish soil layers. However, microbial community composition 
could be affected by different soil hydrothermal and aeration conditions with changes in soil 
layers (Brockett et al., 2012). Due to the harsh environmental conditions in the desert, nutrient 
content and microbial biomass in the deeper soil (e.g., 10—20 cm) can be even higher than those in 
the surface soil and are worth investigating (Rodriguez-Zaragoza et al, 2005; Shamir and 
Steinberger, 2007; Fang et al., 2013). Therefore, microorganisms in the subsurface soil may not 
only respond to precipitation changes but also play a more important ecological role. In some 
hyperarid deserts, such as the Atacama Desert, research involving in soil microbes tends to avoid 
the surface soil with extreme environments and instead targets the deeper soil (Orlando et al., 
2010; Jones et al., 2017; Neilson et al., 2017). However, studies on how the deeper soil microbes 
respond to precipitation and how their responses differ from those of surface soil are rare. 

Alxa Desert in Northwest China seems to be an ideal research area, as it is a transitional zone 
between steppe and typical desert zones, and possesses distinct transitional features (Xiao et al., 
2017). Mean annual precipitation (MAP) declines from east to west, decreasing from 200 mm at 
the foot of the Helan Mountains to less than 50 mm in the desert area (Li et al., 2009). However, 
very few studies have focused on soil microbial changes along precipitation gradients in the 
deserts of this region. Therefore, we selected the study area with a MAP gradient from 150 to 50 
mm in the Alxa Desert and used high-throughput 16S rDNA sequencing to analyze the soil 
bacterial community structure (SBCS), diversity and soil bacterial co-occurrence pattern between 
the surface and subsurface soil layers. We also investigated the environmental factors that might 
influence the structure and diversity of soil bacteria. Our hypotheses are as follows: (1) bacterial 
distribution in the two considered soil layers will both significantly change along the precipitation 
gradient; (2) change in bacterial conditions will be greater in the surface soil layer than in the 
subsurface soil layer; and (3) variations of bacteria along the precipitation gradient will differ 
between the two soil layers, with that in the subsurface layer being more consistent with the 
transitional characteristics of local environment. 


2 Materials and methods 


24 Study area 

Alxa Plateau is a part of Inner Mongolian Plateau, located in the western region of the Inner 
Mongolia Autonomous Region, China. Geographically, it includes the Alxa Desert that 
comprises three major deserts, i.e., Ulan Buh Desert, Tengger Desert and Badain Jaran Desert. 
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Xerophytes constitute the main vegetation (Li et al., 2009) in the Alxa Desert. Affected by 
topography, inland rivers and anthropogenic activities, vegetation coverage varies significantly, 
and sand dunes tend to be fixed in a western to eastern direction. Alxa Desret is characterized by 
a continental climate within an intermediate temperate zone (Xiao et al., 2017). The annual mean 
temperature ranges from 6.0°C to 8.5°C. Precipitation is mainly concentrated between July and 
September, and MAP gradually decreases in a southeastern to northwestern direction. 
Evaporation increases from southeast (2400 mm) to northwest (4700 mm). We divided the area 
into three different subregions based on their precipitation: typical desert (Td), steppe desert (Sd), 
and desertified steppe (Ds), with MAPs of approximately 50, 100 and 150 mm, respectively (Li 
et al., 2009; Pei et al., 2011). 


2.2 Sampling method 


About 3 to 4 sampling sites were selected within each MAP-related subregion in mid to late 
August 2017. No precipitation was observed at least 3 d before sampling occurred at each site. 
We divided the sampled soil into layer A (surface, 0-10 cm) and layer B (subsurface, 10-20 cm) 
(Rodriguez-Zaragoza et al., 2005; Shamir and Steinberger, 2007; Fang et al., 2013). Information 
on sampling sites and environmental conditions is provided in Table 1. No biological soil crusts 
were detected in the sampling sites, while wind-blown sand was found to be widely distributed. 
To ensure consistency in other environmental factors, we collected all the samples from 
semi-fixed sand dunes at the edge of the desert, which closely represent the overall feature 
coverage of the study area (Li et al., 2009). At each sampling site, we selected three Nitraria 
tangutorum Bobr. shrubs with normal growth behaviour for sampling. All soil samples were 
collected below the canopy, approximately 15 cm from the base of shrubs. Around three shrubs, 
we collected soil samples 5 times (approximately 250 g) from each of the two layers and mixed 
15 samples (5 samplesx3 shrubs) from each layer to make a composite sample. A total of 20 
composite samples were collected for all 3 subregions (10 sampling sitesx2 soil layers), and then 
these samples were sealed in the sterilized bags and kept at 4°C in the laboratory. Partial samples 
were stored in a freezer at —80°C until high-throughput sequencing was performed, and the 
remaining parts of the soil samples were used to determine environmental parameters after air 
drying. At the time of sampling, we also tested the samples for soil moisture content (SMC, 


Table 1 Location and environmental characteristics of sampling sites 


Sampling Latitude and Elevation : . Coverage : 
Type età longitude (m) Dominant species (96) Soil type 
40°55'22"N, 
Typical Tdi 100952'40"F. 1031 Nitraria tangutorum, Haloxylon 7.01 
desert 41?23'48"N ammodendron; Gray-brown 
(MAP=50 Td2 10292126"E 904 Reaumuria songarica, Calligonum 2.95 calcic soil 
mm) A191T'53^N mongolicum 
Td3 104*08/45"E 794 3.18 
40?2121"N, 
Sd1 104°44'19"E 1278 17.93 
o3 Nitraria tangutorum, Zygophyllum 
Steppe 39?28'13"N, 
Sd2 1089 xanthoxylon 15.45 
105?36'16"E : - 
nci 00 Reaumuria songarica, Caragana ks d 
mm) Sd3 40°03'12"N, 1471 korshinskii, 26.87 E 
103°55'08"E Pearl russianthistle i 
39°37'16"N, 
Sd4 103°07'54"E 1232 12.29 
38°16'34"N, Nitraria tangutorum, Caragana 
Ds1 Sie 1394 g , g 35.66 
Desertified 103°56'35"E korshinskii, 
steppe 38°54'33"N, Oxytropis aciphylla, Reaumuria Gray desert 
(MAP=150 pe 105?30'15"E 1345 songarica, Ammopiptanthus S soil 
mm) 37954'02"N, mongolicus, 
Ds3 105°12'31"E 1416 Zygophyllum xanthoxylon 28.44 


Note: MAP, mean annual precipitation. 
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g/g) and investigated the vegetation conditions. In addition, we also sampled mature leaves from 
the shrubs to analyze the stoichiometric characteristics of N. tangutorum. 


2.3 Soil bacterial DNA extraction and library preparation 


This step was carried out by the Novogene Genome Sequencing Company and the 
cetyltrimethylammonium bromide/sodium dodecyl sulfate extraction methods were used. To 
meet the criteria for subsequent data analysis, we performed the extraction of DNA from each 
sample 3 times to represent biological replicates. Each replicate was used in sequencing and 
library preparation individually. DNA amplification involved 341F 
(5'-CCTAYGGGRBGCASCAG-3’) and 806R (5'-GGACTACNNGGGTATCTAAT-3') primer 
pair, which corresponds to the 16S V3—V4 regions in bacteria (Michelsen et al., 2014). The 
following polymerase chain reaction (PCR) protocol was used (30.0 pL): (2x) Phusion® 
High-Fidelity PCR Master Mix (15.0 pL, New England Biolabs, USA), 2 pM F-primers (1.5 pL), 
2 pM R-primers (1.5 pL), 10.0 pL template DNA (1 ng/pL), and ultrapure H20 (2.0 pL). PCR 
was performed as follows: initial denaturation at 98°C for 1 min, followed by denaturation at 
98°C for 10 s, annealing at 50°C for 30 s, extension at 72°C for 30 s (30 repetitions), and 72°C 
for 5 min. Electrophoresis was conducted to separate PCR products in 2% agarose gel, and target 
DNA was recovered by cutting strips from the gel. PCR products were purified using the 
GeneJETTM Gel Extraction Kit (Thermo Fisher Scientific, USA). 

The construction of library was carried out by using the Ion Plus Fragment Library Kit 
(Thermo Fisher Scientific, USA). The library was established following quantification and 
testing by applying the Qubit® 2.0 Fluorometer (Thermo Fisher Scientific, USA), and the Ion 
S5TM XL System (Thermo Fisher Scientific, USA) was used for sequencing, generating 
approximately 465 bp single-end reads. We used Cutadapt (v1.9.1) to filter low-quality reads 
(Martin, 2011), and then assigned single-end reads to samples based on their unique barcodes 
and truncated by cutting off the barcode and primer sequences. Splitting sequences were referred 
to as raw reads. To obtain high-quality clean reads, we used the UCHIME algorithm to compare 
the reads with those in the species annotation database to detect and remove chimaeric sequences 
(Edgar et al., 2011; Haas et al., 2011). 


2.4 Determination of environmental factors 


To evaluate the environmental factors, we measured SMC, pH, vegetation coverage (COV), soil 
nutrient factors (contents of soil organic carbon (SOC), total nitrogen (TN), total carbon (TC), 
available nitrogen (AN), total phosphorus (TP), available phosphorus (AP) and available 
potassium (AK)), soil salinity-related factors (contents of the ions SO; , CI, HCO;;, Na*, Mg”, 
Ca?*, K* and electric conductivity (EC)), N. tangutorum stoichiometric factors (total carbon 
(TCn), total nitrogen (TNn), and total phosphorus (TPn)), and soil particle size distribution 
(SPSD). Determination of related environmental factors was carried out in the Key Laboratory of 
Ecohydrology of the Inland River Basin, Chinese Academy of Sciences. 


2.5 Data analysis 


Sequencing data were analysed by Novogene Company (Beijing, China). Sequences with >97% 
similarity were assigned by Uparse (v7.0.1001, http://drive5.com/uparse/) to the same operational 
taxonomic unit (OTU). We used Silva database (v132, https://www.arbsilva.de/) to annotate 
taxonomic information based on Mothur algorithm (Quast et al., 2013). We used Muscle software 
(v3.8.31, http://www.drive5.com/muscle/) for multiple sequence alignment (Edgar et al., 2004) to 
study the phylogenetic relationships between OTUs and dominant species in different samples or 
groups. Finally, we normalized data related to OTU abundance using a sequence quantity 
standard corresponding to the sample with the least sequences and rarefaction. We performed 
subsequent alpha and beta diversity analyses based on the normalized data output. 

Four common alpha diversity indices were calculated with QIIME software (v1.7.0) and used 
to investigate the diversity within a community: Chaol, Shannon (H'), Simpson (D) and 
phylogenetic diversity (PD) whole tree (genetic relationship within a community). For beta 
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diversity, we calculated unweighted UniFrac distance values using QIIME software (v1.7.0) 
based on OTU relationships. Then we used OTU abundance data to further calculate weighted 
UniFrac distance values (Lozupone et al., 2011). In this study, we selected weighted UniFrac 
distance values as the basis for relevant beta diversity analysis. To conduct beta diversity 
analysis, we used analysis of molecular variance (AMOVA, calculated with Mothur software 
(v1.29.0)) to reflect the intragroup differences and assess the significance differences between 
groups (Excoffier et al, 1992; Rivas et al, 2013). Linear discriminant analysis effect size 
(LEfSe, http://huttenhower.sph.harvard.edu/galaxy) was used to identify biomarkers (Segata et 
al. 2011) among MAP-related subregions and soil layers. Default linear discriminant analysis 
(LDA) threshold was set at 24. The Psych and Vegan packages of R software (v2.15.3) were 
used to screen and ordinate environmental factors. We used network topology analysis to reflect 
the differences in soil bacterial co-occurrence patterns among MAP-related subregions and soil 
layers. To remove the influence of community richness differences, we made the numbers of 
bacterial genera in each group consistent (250 genera). The threshold for correlation was 
determined at r>0.80 and P<0.05 conditions. 


3 Results 


3.1 Differences in SBCS and alpha diversity along the precipitation gradient between soil 
layers 


A total of 3,503,236 raw reads were obtained for 60 samples, and 2,949,861 clear reads 
(86.4696—96.0796) were obtained after quality control. After filtering out unknown results, we 
detected 7301 OTUS across 60 samples according to the annotation results. A total of 496 genera 
were annotated, belonging to 201 families, 101 orders, 36 classes and 38 phyla (without 
undefined taxa). The dominant phyla in the soils were Actinobacteria (35.6696 (+13.18%)), 
Proteobacteria (30.5796 (+8.00%)), and Bacteroidetes (13.1696 (+10.79%)) (Fig. 1). The dominant 
genera in the soils were Arthrobacter (5.89%), Sphingomonas (3.78%) and Nocardioides 
(3.32%). 
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Fig. 1 Relative abundance of soil bacterial at phylum level in different MAP (mean annual precipitation)-related 
subregions and soil layers. Td, typical desert (50 mm precipitation); Sd, steppe desert (100 mm precipitation); Ds, 
desertified steppe (150 mm precipitation); A, 0—10 cm soil layer; B, 10—20 cm soil layer. 


Figure 2 shows the statistical results for four commonly used alpha diversity indices. Layer B 
had significantly higher (Wilcoxon rank-sum test, P«0.05) values for all indices than layer A. 
This result indicated that layer B had more species (richness), greater evenness and greater 
evolutionary distance within the community. More importantly, there were obvious differences in 
the changes between the two soil layers along the precipitation gradient. The diversity of layer B 
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increased rapidly and linearly, whereas that of layer A did not have these features. In particular, 
there were no significant differences between Sd.A and Ds.A for all four indices. In terms of the 
differences in SBCS among different treatments (Table 2), the change in SBCS along the 
precipitation gradient between soil layers exhibited some similar characteristics. For example, in 
both layers, difference between Td and Sd was larger than that between Sd and Ds. However, 
unlike the diversity results, the results showed that difference in SBCS along the precipitation 
gradient was greater in layer A than in layer B. 

LEfSe analysis can provide more detailed information on the differences in SBCS among 
subregions in each soil layer (Figs. 3, S1 and S2). Along the precipitation gradient, a total of five 
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1500. E B : 
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Chaol index 
L 
Shannon index 


3 


3 


Simpson index 
z 
& 
PD whole tree 


0.94 


Td Sd Ds , Sd Ds 
MAP-related subregion MAP-related subregion 
Fig. 2 Alpha diversity index in different MAP-related subregions and soil layers. Td, typical desert (50 mm 
precipitation); Sd, steppe desert (100 mm precipitation); Ds, desertified steppe (150 mm precipitation); A, 0-10 
cm soil layer; B, 10-20 cm soil layer; PD, phylogenetic diversity. The lowercase letters indicate significant 
differences among different treatments at P<0.05 level. 


Table 2 Result of analysis of molecular variance (AMOVA) in different MAP-related subregions and soil layers 


Group F value P value Group F value P value 
Ds-Sd 3.021 0.008 Td.A-Sd.A 16.727 «0.001 
Sd-Td 12.343 <0.001 Td.A-Ds.A 15.883 «0.001 
Ds- Id 12.203 «0.001 Sd.A-Ds.A 4.081 0.004 
Td.A-Td.B 5.069 0.007 Td.B-Sd.B 4.615 «0.001 
Sd.A-Sd.B 7.222 «0.001 Td.B-Ds.B 7.022 «0.001 
Ds.A-Ds.B 17.155 «0.001 Sd.B-Ds.B 2.925 0.014 
A-B 13.404 «0.001 


Note: Td, typical desert (50 mm precipitation); Sd, steppe desert (100 mm precipitation); Ds, desertified steppe (150 mm precipitation); 


A, 0—10 cm soil layer; B, 10—20 cm soil layer. 


phyla were identified as the major phyla. Td subregion was characterized by a high abundance of 
Bacteroidetes and Firmicutes. In contrast, Sd and Ds subregions were characterized by a high 
abundance of Actinobacteria, Acidobacteria and Chloroflexi. It is remarkable that, compared with 
layer A, there were only a few biomarkers in layer B in Sd subregion. The abundance of the most 
major differential phylum in layer B was more in line with the change in MAP and showed a 
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continuous increase or decrease along the precipitation gradient. In layer A, however, only 
Firmicutes exhibited this characteristic. Comparing the two soil layers, we found a higher 
abundance of Proteobacteria and Bacteroidetes in layer A than in layer B (Fig. 4). In contrast, 
layer B had a higher abundance of Chloroflexi and Acidobacteria, which was similar to the 
characteristics of Sd or Ds subregion. 
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Fig.3  Cladograms generated according to LEfSe analysis for the precipitation divisions in layer A (a) and layer 
B (b). Biomarkers were statistically significant at LDA 24. Cladograms indicating the phylogenetic distribution of 
bacterial lineages associated with the three subregions. Phylum, class, order, family, genus and species levels are 
listed in order from inside to outside of the cladogram. Red, green and blue circles represent the bacteria enriched 
in Ds, Sd and Td subregions, respectively, whereas yellow circle represents the taxa showing no significant 
difference among the three subregions. Labels for biomarkers for the class, order and family levels are 
abbreviated with a single letter. For the genus and species levels, please see Figures S1 and S2 for details. 


Based on canonical correlation analysis (CCA) and after the screening of environmental factors 
according to VIF analysis, we identified 9 environmental factors (Table 3) showing the highest 
correlation with SBCS in each soil layer (Fig. 5). Results showed that the effects of major 
environmental factors were similar between two soil layers, and the main environmental factors 
can be divided into two distinct clusters. One cluster including MAP, COV and AN was positively 
correlated with the abundance of dominant genera in Sd and Ds subregions (such as Nocardioides 
and Microvirga). The other cluster including salt ions (Ca?* and CI), TP and VCS was positively 
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correlated with the abundance of dominant genera in Td subregion (such as Rhodococcus and 
Pontibacter). However, correlations between environmental factors and SBCS were different 
between two soil layers. Environmental factors that had higher correlations with SBCS in layer B 
were mainly soil salt ions and soil EC, while SBCS in layer A had higher correlations with more 
types of environmental factors, such as MAP, coverage, available AN, VCS and N. tangutorum 
stoichiometric factors (e.g., TCn) (Table 4). 
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Fig. 4 Cladogram generated according to LEfSe analysis for the soil layers and biomarkers showing statistical 
significance at LDA >4. Red and green circles represent the bacteria enriched in the A and B soil layers, 
respectively. 


Furthermore, correlation between environmental factors and diversity differed between two soil 
layers (Fig. S3), and only Ca?* was significantly (P<0.01) correlated with all the diversity indices 
of both soil layers. Interestingly, many environmental factors that were strongly associated with 
SBCS of layer A (e.g., pH, MAP, VCS and COV), were significantly related to the diversity of 
layer B (Table 4; Fig. S3). By contrast, environmental factors that were strongly associated with 
SBCS of layer B (e.g., Ca?', SO? and Cl), were still significantly related to the diversity of the 
same layer. This finding indicated that variation pattern of SBCS and diversity along the 
precipitation gradient was quite different in layer A, while in layer B, it was relatively consistent. 


3.3 Differences in co-occurrence pattern between soil layers 


Based on the relative abundance data at the genus level, we used network topology analysis to 
reflect the differences in soil bacterial co-occurrence patterns among the MAP-related subregions 
and soil layers. The results (Table 5) demonstrated that layer B showed a more obvious and 
monotonous trend than layer A, clustering coefficient (CC) and graph density (GD) increased 
with increasing precipitation, and the average path length (APL) continued to decrease. This 
result indicates that a greater amount of precipitation can support denser and more tightly 
connected communities. Considering that the richness (Chao1 index) increased with increasing 
MAP, this effect may actually be magnified. However, with the increase in precipitation, 
proportion of negative relationship (PNR) between nodes also continued to increase. This finding 
suggests that a greater amount of precipitation also leads to increased intergeneric competition in 
layer B. However, in layer A, changes in CC, GD and PNR were not consistent with the trend of 
MAP. In particular, Sd.A treatment with moderate precipitation had the highest CC and highest 
PNR but the lowest GD (Table 5). This result can be interpreted as the intergeneric competition is 
stronger and the range of close intergeneric associations is strictly limited in layer A in Sd 
subregion. 

In terms of the correlation between environmental factor and intergeneric interaction (Fig. S4), 
we found that MAP and COV were significantly correlated with all network parameters, TP and 
salt ions were significantly (P«0.05) correlated with multiple network parameters. It is worth 
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mentioning that soil moisture likely had no significant effect on SBCS and bacterial diversity in 
this study (Mantel test, r«0.10), but it still showed a strong and significant (P<0.01) correlation 
with the variables describing the internal relationships of soil bacterial community (GD and APL). 
Correlation between GD and soil moisture (r=0.652, P<0.01) was even higher than that between 
GD and MAP (r=0.545, P«0.05). 


Table 3 Differences in environmental factors (soil nutrients, moisture, salinity, pH, particle size distribution, N. 
tangutorum stoichiometric characteristic and vegetation coverage) in different MAP-related subregions and soil 
layers 


Index Td. A Td.B Sd.A Sd.B Ds.A Ds.B 
SOC (g/kg) 1.11+0.71° 0.89+0.35° 1.09+0.53° 0.84+0.21° 1.62+1.15° 1.33+0.86° 
TN (g/kg) 0.16+0.06° 0.21+0.05° 0.26+0.10® 0.2120.10?^ 0.3120.11? 0.20::0.10?^ 
TC (g/kg) 3.85+0.66* 4.76+0.53° 5.87+4.16? 4.10+1.99° 4.61+4.02° 4.07+3.78° 
TP (g/kg) 0.28+0.07° 0.30+0.11° 0.25+0.03? 0.21+0.05° 0.16+0.04° 0.15+0.02° 
AN (mg/kg) 15.10+9.82  16,.28+12.12> 32.39+9.59" 33.95+3.69" 27.54+10.5 25.27+9.11% 
AP (mg/kg) 1.58+0.20° 1.42+0.54° 2.15+0.12* 2.25+1.91® 3.10+1.03" 1.58+0.63° 


AK (mg/kg) 99.60-31.33"  113.43416.60° 155.63-40.83?  112.05-39.01' = 127.27+47.86 135.57+46.21® 


SMC (%) 1.33+0.78° 2.47+0.01° 1.30+0.60° 2.66+1.00* 3.07+0.76* 3.92+1.55° 
SO} (cmol/kg) 4.77+4.75° 10.31+5.83? 1.9543.31° 1.9243.31° 0.12+0.10° 0.09+0.06° 
Cl (cmol/kg) 17.55+28.02° 5.7646.24° 0.50+0.10° 0.15+0.10°° 0.07+0.01° 0.07+0.01° 
HCO; (cmolkg) ^ 0.27+0.06" 0.23+0.04° 0.36+0.08® 0.38+0.08" 0.38+0.06" 0.39+0.06" 
Ca?* (cmol/kg) 2.97+2.09° 7.7144.367 1.90+3.07°° 1.82+3.07"° 0.32+0.08° 0.27+0.05° 
Mg” (cmol/kg) 0.63+0.80" 0.44+0.29" 0.18+0.13° 0.17+0.13° 0.11+0.08° 0.12+0.04> 
K* (cmol/kg) 0.13+0.11? 0.12+0.06° 0.04+0.04° 0.04+0.04° 0.04+0.02° 0.05+0.01° 
Na* (cmol/kg) 15.78+25.01? 6.10+6.29* 0.34+0.21° 0.32+0.21° 0.02+0.03° 0.05+0.02° 
EC (dS/m) 8.04411.31° 5.06:3.86? 1.36+6.14" 4.34+6.15° 0.31+0.03° 0.39+0.05° 
0] pH —— 7443016" 7381005  7.67£0.16°  7.80£0.16° 7.722008" —— 7.702025 — 
| VCS6) 5913271 —— 3224362 0221044: —— 0,033007!  0.7240.66° ——— 0334029: — 
CS (96) 14.30+5.99? 8.89+6.82” 3.84+1.02> 5.95+5.24> 21.07+18.57" 13.05+11.79° 
MS (%) 21.58+9.17° 21.28+8.50° 26.69+8.40° 25.57+15.04° — 27.67+17.71° 27.31+16.48° 
FS (%) 43.32+7.39° 48.24+7.74? 47.26:10.84^ 50.13+21.77? 32.89:18.77  40.77+12.22° 
VFS (96) 11.13+5.53” 12.66+3.23? 8.78+3.91” 7.90+5.14° 11.43+9.37” 11.47+6.58” 
Silt (%) 3.20+1.61° 4,82+2.31° 11.20+11.19* 9.10+9.37* 5.33+7.42# 6.20+7.99* 
Clay (%) 0.57+0.50* 0.88+0.67° 2.032.278 1.32+1.67° 0.89+1.55° 0.87+1.51° 
EE TCn(gkg ^ 3847114689 —  37766:4220 377233759 
TNn (g/kg) 18.3921.92* 20.45 + 7.31? 21.3246.94? 
TPn (g/kg) 1.23+0.08” 1.52 + 0.78? 1.15+0.09° 
DM COV(6) 438498  1814t56^ 38954089 


Note: SOC, soil organic carbon; TN, total nitrogen; TC, total carbon; TP, total phosphorus; AN, available nitrogen; AP, available 
phosphorus; AK, available potassium; SMC, soil moisture content; EC, electric conductivity; VCS, very coarse sand (2.00-1.00 mm); 
CS, coarse sand (1.00—0.50 mm); MS, medium sand (0.50-0.25 mm); FS, fine sand (0.25-0.10 mm); VFS, very fine sand (0.10—0.05 
mm); silt, particle size is 0.050—0.002 mm; clay, particle size is «0.002 mm; TCn, total carbon of N. tangutorum; TNn, total nitrogen of 
N. tangutorum; TPn, total phosphorus of N. tangutorum; COV, coverage. Td, typical desert (50 mm precipitation); Sd, steppe desert 
(100 mm precipitation); Ds, desertified steppe (150 mm precipitation); A, 0—10 cm soil layer; B, 10-20 cm soil layer. The lowercase 
letters are Wilcoxon's rank sum test results, and the lowercase letters within the same row indicate significant differences among 
different treatments at P<0.05 level. Mean+SD. The abbreviations are the same in Figs. 5, $3, S4 and Table 4. 
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Fig. 5 Canonical correlation analysis (CCA) ordination of bacterial community structure at the genus level in 
relation to major environmental factors in layers A (a) and B (b) 


Table 4 Mantel test comparing the correlation between environmental factors and soil bacterial community 
structure in layer A and layer B 


Layer A r value P value Layer B r value P value 
VCS 0.7811 <0.001 Ca” 0.7993 <0.001 
Na* 0.6744 «0.001 SO; 0.7242 <0.001 
pH 0.6366 <0.001 cr 0.6586 <0.001 
MAP 0.6125 <0.001 Na* 0.6384 «0.001 
Cr 0.5879 «0.001 HCO4 0.5660 «0.001 
SOi 0.5877 «0.001 K* 0.4798 «0.001 
EC 0.5766 «0.001 TP 0.4403 «0.001 
Ca” 0.5717 <0.001 Mg?* 0.4013 0.002 
Mg” 0.4498 <0.001 pH 0.3987 <0.001 
K* 0.4090 «0.001 EC 0.3601 «0.001 
AN 0.3987 «0.001 VCS 0.3558 <0.006 
COV 0.3099 0.003 MAP 0.3397 <0.001 
TCn 0.2961 <0.001 AN 0.2644 <0.001 
AP 0.2855 0.004 COV 0.2077 0.009 
VFS 0.2020 0.008 


Table 5 Parameters of network topology analysis (genus level) in different MAP-related subregions and soil 


layers 
Type MD CC GD APL PNR (%) 
Td.A 0.6502 0.4980 0.0242 4.443 24.41 
Sd.A 0.5014 0.5539 0.0230 4.272 39.15 
Ds.A 0.3332 0.5299 0.0354 3.691 25.14 
^ —  TàB 05215 —— 0520 | 0.0332 | A072. o 28.19 
Sd.B 0.5512 0.5407 0.0489 3.576 36.01 
Ds.B 0.3704 0.5721 0.0527 3.303 40.14 


Note: Community richness of each group had been modified to be consistent (250 genera). MD, modularity; CC, clustering coefficient; 
GD, graph density; APL, average path length; PNR, proportion of negative relationship. The abbreviations are the same in Figure S4. 
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4 Discussion 


Different bacterial variation characteristics between soil layers along a precipitation gradient can 
represent the effects of soil depth on the soil microbial response to climate change. Such an effect 
may change the degree and trend of the bacterial response to precipitation change. In this study, 
there were obvious differences between soil layers with the changes in soil bacterial conditions 
along the precipitation gradient. The results showed that in layer B, with the increase in 
precipitation, soil bacterial diversity, connectivity within community and intergeneric negative 
interactions also increased (Fig. 2; Table 5). The change in SBCS was also more consistent with 
the change in diversity in that layer. In contrast, SBCS of layer A changed more dramatically 
along the precipitation gradient than that of layer B (Table 2), but soil bacterial diversity changes 
were limited (Fig. 2). Furthermore, changes in most parameters of co-occurrence patterns in layer 
A did not fit the changes in MAP. The above results are consistent with our first and third 
hypotheses but partly not with the second hypothesis. The differences in environmental factors 
between the two soil layers were mainly nonsignificant (Table 3), and the effects of the main 
environmental factors on SBCS were also similar (Fig. 5). Nevertheless, soil moisture may be an 
important factor that causes different soil bacterial variation characteristics along the precipitation 
gradient between the two soil layers. 

In general, differences in precipitation affect soil microorganisms by influencing changes in 
soil moisture (Zhang et al., 2018). However, correlation between soil moisture and 
SBCS/diversity in our study was very low (r«0.10). This result should be due to the differences in 
soil salinity (Table 3), which distorted the biological availability of soil moisture (Zahran, 1997; 
Rath et al, 2017) among subregions. Even so, we found that SMC was still the only 
environmental factor with significant differences between soil layers throughout the study area, 
and the biomarkers in layer B were similar to those in the humid subregions (Figs. 3 and 4). 
Therefore, soil moisture may still be one of the main factors leading to the difference in soil 
bacterial conditions between the two soil layers. More importantly, we found that there was still a 
high correlation between SMC and connectivity within bacterial community. This finding should 
be related to the strong control of soil moisture over the diffusion of microbial resources 
(Banerjee et al., 2016). Our results show the particularity of the co-occurrence patterns of soil 
bacteria when faced with changes in SMC and show the different effects of MAP and SMC on 
soil bacteria. This result helps to explain the reason for the different characteristics of soil bacteria 
between different soil layers along the precipitation gradient. 

Due to the high heterogeneity of abiotic environments in soil, there are many disconnected 
pores that filled by air or water (Kuzyakov and Blagodatskaya, 2015). As soil becomes drier, 
there is less water in soil pores, resulting in disconnected resource islands (Schimmel, 2018). For 
drier desert surface, there is less water in soil pores and may create more disconnected soil niches. 
This condition may affect the community structure, diversity and co-occurrence patterns of soil 
bacteria at the same time. On the one hand, such relatively isolated niches may provide shelter for 
different types of soil microorganisms and lead to increased diversity (Šťovíček et al., 2017). 
However, selection pressures generated by extreme surface environmental conditions (e.g., large 
temperature differences and strong ultraviolet radiation) tend to reduce the diversity of layer A, 
and these environmental conditions do not improve with increasing precipitation within a certain 
range. Due to the combined influence of the above factors, diversity of soil bacteria in layer A did 
not change significantly with precipitation in this study. On the other hand, the disconnected 
niche formed by a low SMC can weaken the intercommunity connectivity and may further 
weaken the stability of SBCS in layer A. As a result, SBCS in layer A changed more dramatically 
than that in layer B with the changes in precipitation. 

Vegetation is closely related to soil microorganisms. Although the vegetation types were 
consistent in this study, vegetation can still indirectly affect soil microorganisms by changing the 
soil environment. For example, shrub mulching can increase soil infiltration capacity and reduce 
evaporation in shallow soil (Scholes et al., 1997). It has been suggested that such retention of soil 
moisture by vegetation may be most significant in areas with a low precipitation (e.g., Td 
subregion) (D'Odorico et al., 2007). This condition can explain why the moisture content in Sd 
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subregion was slightly lower than that in Td subregion in layer A (Table 3) and may further result 
in the more unique soil bacterial co-occurrence pattern observed in Sd.A treatment (Table 5). In 
addition, vegetation itself also responds to changes in precipitation and changes the impact on soil 
microorganisms. In arid regions, there is a strong competition between soil microorganisms and 
aboveground vegetation for nitrogen and phosphorus due to limited soil nutrients (Cui et al., 
2018). With increasing precipitation, the competitiveness of vegetation for nutrients can be 
enhanced (Dijkstra et al., 2015). In this study, soil available nitrogen and total phosphorus were 
the few soil nutrient factors with a high correlation with SBCS (Table 4). Contents of soil AN and 
TP in Ds subregion were even lower than those in Sd subregion with increasing vegetation 
coverage (Table 3). Such nutrient limitations may offset the environmental gain from increased 
precipitation on soil bacteria and lead to a small difference in SBCS between Sd and Ds 
subregions. Notably, Alxa Desert is one of the main invasion channels of the East Asian winter 
monsoon (Yao et al., 2011), and most litter can be blown away or buried as a result of intense 
blown sand activity. For other deserts, the accumulation of litter under the canopy or the presence 
of biological crust may result in a significant difference in the extent of such nutrient limitation 
between layers. 

For long-term climate change, the stability of soil environmental conditions will also affect the 
response degree of soil bacteria. A more stable soil environment will allow more time for soil 
bacteria to acclimate. For the surface soil, intense evaporation and rapid infiltration lead to a rapid 
SMC change. Such a rapid change in SMC, or called re-wet, can lead to a short period of 
increasing mortality and decreasing diversity of soil bacteria (Schimmel, 2018). For the deeper 
soil, the drying of the upper layers prevents soil moisture migration and further evaporation 
(Noy-Meir, 1973; Epstein et al., 1997). Therefore, the deeper soil can maintain a higher SMC and 
more stable moisture condition than surface soil. However, in addition to increasing drought, 
climate change can lead to more severe and erratic precipitation events (Easterling et al., 2000). 
Severe precipitation events are beneficial to supplementing deep SMCs and maintaining them for 
a relatively long time. Considering that desert soil microorganisms have evolved various 
physiological strategies to resist drought (Barnard et al., 2013), short-term SMC improvements 
can also lead to a long-term response of soil bacterial community. Thus, how changes in 
precipitation patterns affect the stability of soil environment at different depths and how microbes 
respond at different depths still need further study. 


5 Conclusions 


In this study, soil bacterial status at both the surface and subsurface was found to significantly 
vary along the precipitation gradient. Differences in soil salinity and nutrient limitation affected 
the difference degree in SBCS in adjacent subregions and can affect both soil layers. However, 
variation in each soil layer along the gradient also differed: response of SBCS was more intense 
in layer A, while layer B was more consistent with the precipitation gradient in terms of diversity 
and co-occurrence patterns. Layer B can support communities with higher diversity and closer 
internal relationships but more internal competition than layer A. The changes in these variables 
were more in line with the changes in MAP. Therefore, it is necessary to distinguish the response 
of surface and subsurface soil bacteria to changes in precipitation. In addition, soil moisture 
content may play multiple roles in mediating environmental conditions and soil bacterial 
community characteristics. When predicting the trends in desert soil bacterial conditions 
associated with precipitation, it is necessary to distinguish surface and subsurface soils. The 
trends of bacterial condition in surface soil needs to consider the influence from more types of 
environmental factors, and the specific depths at which these factors can affect are need to 
specific research. 
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Appendix 
EE Ds.A ESA — NEMTUA 


p_Bacteroidetes 
c_Bacteroidia 
o_Cytophagales 
p_Firmicutes 

g Pontibacter 

c Bacilli 

c Gammaproteobacteria 
o Bacillales 

f Planococcaceae 

f Hymenobacteraceae 

f Cyclobacteriaceae 

o Flavobacteriales 

f Flavobacteriaceae 

f Rhodobacteraceae 

o Rhodobacterales 

s Rhodococcus erythropolis 
g Rhodococcus 

f Nocardiaceae 

g Salinimicrobium 

c Longimicrobia 

f Longimicrobiaceae 

o Longimicrobiales 

p. Actinobacteria 

c unidentified Actinobacteria 
p. Chloroflexi 

f. Chitinophagaceae 

o Chitinophagales 

p. Acidobacteria 

f Nocardioidaceae 

o unidentified Gammaproteobacteria 
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Fig. S1 Histogram of linear discriminant analysis (LDA) score distribution shows the categories with LDA 
score greater than the set value (7 4) in layer A (0—10 cm) along the precipitation gradient. The length of the 
histogram indicates the magnitude of the influence of biomarkers. 
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Fig. S2 Histogram of linear discriminant analysis (LDA) score distribution shows the categories with LDA 
score greater than the set value (> 4) in layer B (10-20 cm) along the precipitation gradient. The length of the 
histogram indicates the magnitude of the influence of biomarkers. 
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Fig. S3 Major environmental factors associated with alpha diversity in layers A (a) and B (b). Environmental 
factors are significantly correlated with at least one alpha diversity index. Asterisks indicate significance of the 
Spearman's r correlation. *, P«0.05 level; **, P«0.01 level. 
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Fig. S4 Major environmental factors associated with parameters of network topology analysis. Environmental 
factors are significantly correlated with at least one parameter. Asterisks indicate significance of the Spearman's r 
correlation. *, P<0.05 level; **, P<0.01 level. 


